arXiv: 1509.02397v2 [astro-ph.GA] 6 Mar 2016 


Mon. Not. R. Astron. Soc. 000 , 000-000 (0000) Printed 8 March 2016 (MN lAT^^X style file v2.2) 


Stellar Dynamics around a Massive Black Hole I: 
Secular Collisionless Theory 

S. Sridhar^’^ and Jihad R. Touma^’^ 

^ Raman Research Institute, Sadashivanagar, Bangalore 560 080, India 

^ Department of Physics, American University of Beirut, PO Box 11-0236, Riad El-Solh, Beirut 11097 2020, Lebanon 
® ssridhar@rri.res.in ^ jtOO&aub.edu.lb 


8 March 2016 


ABSTRACT 

We present a theory in three parts, of the secular dynamics of a (Keplerian) 
stellar system of mass M orbiting a black hole of mass M, ^ M. Here we 
describe the collisionless dynamics; Papers II and III are on the (collisional) 
theory of Resonant Relaxation. The mass ratio, e = M/M, <C 1, is a natu¬ 
ral small parameter implying a separation of time scales between the short 
Kepler orbital periods and the longer orbital processional periods. The colli¬ 
sionless Boltzmann equation (CBE) for the stellar distribution function (DF) 
is averaged over the fast Kepler orbital phase using the method of multiple 
scales. The orbit-averaged system is described by a secular DF, F, in a re¬ 
duced phase space. F obeys a secular CBF that includes stellar self-gravity, 
general relativistic corrections up to 1.5 post-Newtonian order, and external 
sources varying over secular times. Secular dynamics, even with general time 
dependence, conserves the semi-major axis of every star. This additional in¬ 
tegral of motion promotes extra regularity of the stellar orbits, and enables 
the construction of equilibria, Fq, through a secular Jeans theorem. A lin¬ 
earized secular CBE determines the response and stability of Fq. Spherical, 
non-rotating equilibria may support long-lived, warp-like distortions. We also 
prove that an axisymmetric, zero-thickness, flat disc is secularly stable to all 
in-plane perturbations, when its DF, Fq, is a monotonic function of the an¬ 
gular momentum at fixed energy. 

Key words: galaxies: kinematics and dynamics — galaxies: nuclei — Galaxy: 
center 


1 INTRODUCTION 


Massive black holes (MBHs), consuming gas and s tars, are beheved to be r esponsible for the radiant 
power of quasars and other active galactic nuclei (|Reeslll984l : lKroliklll9991 1 . There is also dynamical 
evidence for MBHs with masse s M, in the range 10^ — 10^^ Mq in the nuclei of more nearby 
galaxies ( Kormendv fc Ho 201,11 ) . Coevolution (or otherwise) of a MBH and its host galaxy has 
been discussed using rnany s tudies of correlations of M, with galaxy properties (IKormendv &: Ho 


2013; 


Heckman fc Best 2014l i . Ground-based observations provided the first dynamical evidence for 
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2 Sridhar & Touma 


MBHs in nearby galaxies (jKormendv fc Richstonelll995l ). They also revealed nuclear star clusters 
where stellar number densities far exceeded those of globular clusters. Early work on the stellar 
dynamics of galactic nuclei considered several aspects of the interactions bet ween the MBH and its 
stellar environment: fee ding of stars to a MBH through two-body relaxation (jBahcall fc WolfUlQTfil : 
Cohn fc KulsrudI 1 197811: th e effect of the adiaba tic growth of a MBH on the formation of stellar 
density cusps ( YoungHigSol l. and on stellar orbits ( Goodman &: Binnevlll984); destruction by a MBH 
of the regularity of stellar orbits in the inner parts of triaxial galaxy ( Gerhard &: Binnev IQSsI l. 
With the Hubble Space Telescope spatially resolved kinematics, of the dense nuclear clusters of stars 
populating the MBH’s sphere of influence, became available. Recent ground-based observations 
with adaptive optics have resolved the sphere of influence of the MB H at our Galactic Gentre so well, 


that the orbits of many individual stars are now being followed (jGenzel. Eisenhauer fc Gillessen 

2ninl l. 


Our interest is in these dens e stellar s y stems in galactic nuclei lying within the sphere of 
influence of the MBH. Eollowing Tremaine ( 20051 ) we may consider a stellar system orbiting a 
MBH to be Kepleriai^ when the gravitational force, over a large range of radii, is dominated by 
the MBH, and the dominant non-Keplerian force arises from the self-gravity of all the stars. In 
the absence of self-gravity the orbit of a star follows a closed Kepler ellipse around the MBH, 
because the radial and azimuthal frequencies are equal to each other (i.e. they are “degenerate” 
or “resonant”). When the small self-gravity of the cluster is taken into account, the orbit of a 
star is usefully represented by its instantaneous osculating Kepler ellipse, which evolves slowly over 
secular time scales that are much longer than the fast orbital times around the MBH. Progress in 
the under s tandi ng of the dynamics of Keplerian stellar systems can perhaps be traced to two papers: 
Tremainel (1199511 on the st ructure of the (then) double nucleus of the Andromeda galaxy (M31), and 
Rauch fc Tremaine ( 19961 ) on the role of resonant relaxation in stellar systems with degenerate (or 
resonant) frequencies. The former has stimulated much work on (collisionless) dynamics whereas 
the latter proposed a new and efficient (collisional) relaxation process called Resonant Relaxation, 
which is the statistical mechanics of Keplerian stellar systems. Two decades of research has paved 
the way for the construction of triaxial, even lopsided, equilibria around MBHs, and simplified 
the modeling, stability analysis, and numerical simulation of nuclear star clusters. We now have 
a better appreciation of the dynamics and the statistical mechanics of Keplerian star clusters, 
and of the additional role of the MBH as a consumer of stars. But the whole subject lacks a 
clearly formulated and general mathematical foundation. Our aim is to fill this gap and provide 
a secular theory of the structure and evolution of Keplerian stellar systems, described by the slow 
time evolution of a secular distribution function (DE), defined in a reduced phase space. Here 
(Paper I) we derive the (collisionless) theory descr ibing their dynamics, begin ning with the standard 
description of co llisionless stellar systems given in Binnev &: Tremaine ( 2008l l. In Sridhar &: Touma 
(12015: Paper Ill ) we deal with their (coll isional) s tatist ical mechanics and presen t a first-principles 


theory of resona nt relaxation based on Gilbert ( lOGsl ). This is then applied in 


Sridhar fc Touma 


(12016: Paper Hil l to the physical k inetics of axisym metric, zero-thickness, flat discs. 

The sketch of a model offered in lTremaind (jl995l l — for a lopsided Keplerian disc of stars around 
a MBH in the nucleus of M31 — has since been fleshed out in a number of numerical models (Statlei 


1999; 

Bacon et al. 

2001; 

Salow &: Statler 

2001: 

Sambhus & Sridhar 

2002; 

Peiris Sz Tremainel 

2003; 

Salow k Statler 

2004; 

Kazandiian & Touma 

2013; 

Brown & Mae:orrian 

2013). This has brought to 


the forefront deep and interesting questions regarding the general nature of secular collisionless 
evolution. Below we give a brief time-line of developments in the secular collisionless dynamics of 


1 


Tremaine J2005h uses the term “near-Keplerian”; 


we have shortened it to Keplerian. 
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Stellar Dynamics around a Massive Black Hole I 3 


Keplerian stellar systems. Work in this field may be broadly divided into explorations of orbital 
dynamics, construction of stationary DFs, linear stability of stationary DFs, and the nonlinear 
evolution of the DF. 

Orbital Dyna mics: The orbital s tructur e o f model galactic-type pot e ntials with a centra l MBH 
was studied in Sridhar fc Touma ( 1997al l]): Merritt fc Valluri 1 19991 ). Sridhar fc Touma ( 19991 ) 
imported the Gaussian orbit-averaging methods of celestial mechanics into stellar dynamics within 
the sphere of influence of a MBH. In this procedure a stellar orbit is averaged over the (fast- 
varying) orbital phase of the corresponding osculating Kepler ellipse, a direct consequence of which 
is the near conservation of the semi-major axis of the stellar orbit. Since there arises an extra 
(approximate) integral of motion, they concluded that a MBH promotes extra order on stellar 
orbits within their sphere of influence. Orbit-averaging takes us to a reduced phase space describing 
secular dynamics: this is the theatre of the coupled evolution of the orbital angular momentum 
and Lenz vectors — giving the eccentricity and spatial orientation of the orbit — keeping fixed the 
semi-major axis. Studying orbits in planar, lopsided model disc potentials, they found a family of 
stable orbits that were lopsided in the same sense as the disc potential. The more complex orbital 


struct u re of triaxial syst e ms w as also studied in a number of papers: see e.g. ISambhus &: Sridhar 


( 2nnd i: Poon fc Merritt! ( 20011 ): Merritt fc Vasiliev ( 2011 ). Of particular i nterest is a fa mily of 
centrophilic orbits that bring stars close to the MBH. These are reviewed in Merritt ( 201,‘ll ). which 
also offers a treatment of the post-Newtonian (relativistic) effects of the MBH on stellar orbital 
dynamics. 

Linear Stability of Stationary DFs: Stationary DFs are often taken as first approximations 
for describing the nuclear star clusters, and their stability to small perturbations is an important 
question. This is, in general, a difficult problem because the long-range nature of gravitational 
interactions leads to an integral eigenvalue problem. Two types of stellar systems have been studied: 
(a) Zero-thickness fiat discs, and (b) Spherical non-rotating systems . 

(a) Linear Stability of Discs: Sridhar. Sver fc Touma ( 19991 ) used the Laplace-Lagrange 
secular theory to study slow, lo psided m = 1 modes (w here m is the azimuthal wavenumber) of 
cold stellar discs around MBHs. Lee &: Goodman ( 1999l i considered nonlinear m = 1 spiral density 
waves in fluid discs, and used a variational metho d to derive the disp e rsion relation and angu¬ 
lar momentum ffuxes. The nu merical sim u lation s of Jacobs fc Sellwo^^ ( 20011 ) provided evidence 
for long-lived lopsided modes. iTremainel (120011 ) demonstrated that a variety (fluid, stellar and 
softened-gravity) of rotating discs were neutrally stable to m = 1 modes in the Wentzel-Kramers- 
Brillouin (WKB) limit. He formulated a linear integral eigenvalue problem for softened-gravity 
discs, and demonstrated that all m = 1 modes were neutrally stabl e; the problem was also so lved 
numerically to determine the eigenfrequencies and eigenfunctions. iJalali fc Tremainel (120121 ) de¬ 
rived the WKB dispersion relationship for slow mo des in a collisi onless disc with a Schwarzschild 
DF and showed that modes of all m were stable. Touma ( 2 OO 2 I 1 presented a Laplace -Lagrange 
theory for discs composed of softened Gaussian wires, confirmed the stability result of iTremaine 
(j200ll ) for prograde discs, and also demonstrated that counter-rotating discs were unstable to slow, 
lopsided m odes: a possible orig i n of the lopsided nuclear disc of M31 may be through such an 
instability. ISambhus fc Sridhar! pOOJ) built on this to include a retrograde population of stars in 
their planar model of the nuclear disc of M31 , and suggested that this population could arise nat¬ 
urally through an accretion event. Tremain^ ( 2005l l considered stellar discs with DFs even in the 
angular momentum, an d empty lo s s con es (i.e. DF is zero at zero angular momentum). Using the 
variational princi ple of iGoodmanI (1198811 he argued that r nany of them are likely to be unstable 
to m = 1 modes. Polyachenko, Polyachenko &: Shukhman ( 2007l l considered mono-energetic discs 
dominated by nearly radial orbits, and found that they are prone to loss cone instabilities of all m. 
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4 Sridhar & Touma 


if there is some amount of counter-rotating stars. Sridhar fc Saini ( 20ld i derived the instability of 
m = 1 modes of softened-gravity, counter-rotating discs in the WKB limit. iGulati. Saini fc Sridhar 
( 2 OI 2 I 1 studied the linear integral eigenvalue problem for these discs and determined eigenvalues, 
precession frequencies and growth rates. 


(b) Linear Stability of Spherical Non rotating Systems: iTremaind ()2005l l formulated 
the eigenvalue problem for the linear secular modes of spherical non-rotating systems and deter¬ 
mined the following: (i) Modes with arbitrary I (where I is the latitudinal wavenumber) are of 
two types — either pairs of modes that are damped and growing in time, or pairs of oscillatory 
modes; (ii) Lopsided I = 1 modes are stable when the DF is either a decreasing function of angu¬ 
lar momentum (at fixed energy), or the DF is an increasing function of angular momentum and 
has a non-empty loss cone; (hi) DFs that are increasing functions of angular momentum with 
an empty l oss cone are neutrally stable to an / = 1 mo de that corresponds to a uniform dis¬ 
placement. Polvachenko. Polvachenko fc Shukhman ( 200?! ) considered a loss-cone instability for 
mono-energetic spherical stellar systems dominated by nearly radial orbits. They found that DFs 
which are non monotonic functions of the angular momentum may be unstable to modes with 
I ^ 3 linear perturbations; a necessary condition for instability is the retrograde precession of 
orbits, which obtains naturally in such clusters. 

Nonlinear evolution of the DF: Numerical studies of the collisionless dynamics of Keple- 
rian stellar systems have been more limited than their galactic counterparts. This is partly be¬ 
cause observations are not yet precise enough to finely resolve the sphere of influence of central 
black holes (the Galactic center and M31’s nucleus being the exceptions), and partly because 
the kind of N-body simulations required to resolve the collisionless regime are computationally 


struct non-axisvmmetric equilibria, both kinematic ( 

Peiris & Tremaine 

20031 and self-consistent 

(Poon & Merritt 

2001 

Sambhus & SridhaJ2002 

Brown &: Magorrian 1201311. Earlv numerical sim- 

Illations ( 

Bacon et al. 

2001; 

Jacobs Sz SellwoodI 

20011 

showed collisionless relaxation of perturbed 
built a numerical algorithm which is partic- 

stellar discs. 

Touma. Tremaine Kr, Kazandiian (20091 


ularly suited to the numerical exploration of secular clusters. This is based on a generalization of 
an orbit-averaging algorithm, first introduced by Gauss, to systems with Plummer-softened inter¬ 
actions. The simulations proceed over secular time steps, in a vectorial framework which avoids 
singularities in the usual Keplerian elements. Averaging of force evaluations, which is costly, is 
ideal for parallel computations on a Beowulf cluster with or without GPU option. Their work 
explored various violent regimes of secular interactions, and provides the first evidence for colli¬ 
sionless relaxation of secularly unstable counter-rotating clusters into lopsided configurations. The 
(fully nonlinear) collisio nless Boltzmann equation (CBE) for secular evolution was presented in 
Touma fc Sridha r! (12012 11. and used t o exp lore the nonlinear, collisionless dynamics of counter- 
rotating discs. iKazandiian fc Touma ( 20131 ) performed N-body simulations of unstable counter¬ 
rotating 3-dim clusters, and demonstrated that they relax collisionlessly into triaxial lopsided 
configurations; these appear to be promising models of the lopsided stellar nucleus of M31. 

We are now in a position to bring all the work together under a general framework, describing 
the secular collisionless evolution of Keplerian stellar systems. We begin in § 2 by casting the 
standard collisionless Boltzmann equation for the phase space DF in the Delaunay action-angle 
variables (which are natural phase space coordinates for the pure Kepler problem). In § 3 the 
method of multiple scales is used to derive an explicit functional form for the DF. Complementing 
this is Appendix A on canonical perturbation theory applied to an individual stellar orbit. In § 4 
self-gravity is augmented by relativistic corrections and external gravitational perturbations, to 
obtain the secular CBE governing the evolution of the secular DF. A secular Jeans theorem is 
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Stellar Dynamics around a Massive Black Hole I 5 


proved in § 5, and used to discuss secular equilibria. Linear secular perturbations and stability are 
discussed in § 6, where some results are presented for two types of systems: (a) Spherical non¬ 
rotating systems; (b) Zero-thickness flat axisymmetric discs. Concluding remarks are offered in 
^ 7. 


2 COLLISIONLESS EVOLUTION OF KEPLERIAN STELLAR SYSTEMS 

In a Keplerian stellar system the gravitational force over a large range of radii is dominated by 
the MBH, and the dominant non-Keplerian force arises from the self-gravity of all the stars. The 
MBH also plays an additional role in destroying stars on orbits with near-zero angular momentum, 
through tidal disruption or being swallowed whole. Let us consider such a system composed of a 
large number (N 3> 1) stars, of total mass M and size R, orbiting a MBH of mass M,. Then we 
must have (1) e = (M/M,) <C 1, and (2) r, = 2GM,/c^ such that {r,/R) <C e; this ensures 
that, over most of the cluster, stellar self-gravity dominates relativistic corrections to the MBH’s 
gravity. Then e is the natural small parameter for the dynamics of non-relativistic star clusters. 

The important time scales in the life of a Keplerian stellar system are: (i) The short Kepler 
orbital timescale, T^ep = 27r ; (ii) The longer secular timescale, Tgec , over which the 

stars behave as if they were Keplerian Ellipses being torqued collisionlessly by the mean gravi¬ 
tational field of the cluster; (iii) The resonant relaxation timescale, Tres, over which the angular 
momenta (but not the energies) of stellar orbits diffuse due to gravitational collisions between stars; 
(iv) The 2-body relaxation timescale, Tj-eiaxi over which both the energies and angular momenta 
of the stellar orbits diffuse, due to gravitational collisions between stars. The different time scales 
are related to each other by, Tkep : T/ec : Res : Relax = 1 : : Ns~^ : Ne~‘^ H For galactic nuclei 

of interest Relax is usually greater than the Hubble time, so the usual 2-body relaxation process 
is not important. However, Res can be less than the Hubble time, so there are three regimes of 
interest: 


• Short times {t < few R;ep Rec ) during which there are variations in the stellar distribution 
over Kepler orbital times. These have indeed been directly observed in our Galactic Centre. The 
general theory is undeveloped and could be important for understanding the first phase of the 
collisionless relaxation of a star cluster, were it to form in the sphere of influence of the MBH. This 
regime is not studied in this paper. 

• Intermediate times (Rep ^ t ~ few Rec ^ Res) over which the cluster behaves collision¬ 
lessly. This secular regime is the focus of the present work (Paper I). 

• Long times (Rec ^ t ~ few Res Relax) over which the system evolves collisionally due to 
resonant relaxation. This is explored in Papers H and HI. 


The two Keplerian systems that have been studied best are also th e closest ones, in which the sphere 
of influence of the MBH is well-resolve d (Kormem y fc 31 boKlI ). The Galactic Centre MBH has 
an estimated mass of ~ 4.6 x 10® Mq ( Yelda et al.ll2014l L and is surrounded by many young and 
old stars, with e ~ 10“^ — 10“^ depending on wha t substructure is b eing studied. The ‘triple 
nucleus’ of M31 has a MBH of mass ~ 1.4 x 10® Mq ( Bender et al.ll2005l ). with a similar values of 
e, depending on where the disc edge is located. One may ask if these values of e are small enough. 
For M31 the answer seems to be in the affirmative because, as discussed earlier, the Keplerian 
disc model provides a natural explanation of the lopsidedness in terms of the trapping of orbits 


^ This estimate of Tres applies to the diffusion of the magnitude of angular momentum, or scalar resonant relaxation. Loga¬ 
rithmic corrections have been neglected in the estimate for Treiax- 
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6 Sridhar & Touma 


by self-gravity. The orbits of many individual stars close to the Galactic Centre MBH have been 
determined, a nd for these t he dy namics is certainly Keplerian. Disc-like structures have also been 
observed, e.g. Yelda et al. ( 2014l i. but secular theory has not progressed to the stage of making 


clear predictions that can be tested. 


2.1 Collisionless Boltzmann eqnation in the Delaunay variables 

Let r be the position of a test star relative to the MBH, t be the time variable, and u = (dr/dt) 
be the relative velocity. A large N stellar system can be described at time t by f{r,u,t), the 
distribution function (DF) in the 6-dim single-particle phase space {r, u}. The DF is positive and, 
when the MBH is not considered a sink of stars, the total mass of the stellar system is conserved. 
Then the DF can be normalized for all time: 

J f{r,u,t)drdu = 1, (1) 


making / a (single-particle) probability distribution function. The gravitational force experienced 
by a star is due to the MBH and all the other stars. Per unit mass the former is —, 
and the latter can be written as —d(e(p)/dr where £ip{r,t) is the mean-field self-gravitational 
potential given by, 

ip{r,t) = -GM, [ 4^:^l^dr'dn'. (2) 

J T — 1" \ 

We write the acceleration of the MBH due to all the stars as £A,{t), where 


A,{t) 


GM, J f{r,u,t)^ 


dr dn. 


( 3 ) 


Note that (p and A, have been defined such that they are 0(1) quantities. Then the relative 
acceleration of a star (with respect to the MBH) is obtained by subtracting the acceleration of the 
MBH from the net gravitational force (per unit mass) on the star: 


du GM, ^ dip 
^ = - ~ ~ 

The stellar system may be regarded as collisionless over times shorter than the relaxation times 
for exchanging energy or angular momentum through stellar gravitational encounters. Then the 
time evolution of the DF is governed by the collisionless Boltzmann equation (CBE): 


d/ _ df dr df du df 
dt dt dt 9r dt du 


dt 


+ u 


dr 


+ 


GM ,. dp 

- V — E - 

r2 Qj, 



du 


= 0 . 


It can be verified that the CBE preserves the normalization of eqn.([T]). Using Poisson Brackets 
(PBs), the CBE can be rewritten compactly as: 

df 

+ [/i -f^org](g) =0, (5) 

where 

, , GM, , , ^ , 

Hovg{r,u,t) = — -h £p{r,t) + £r-A,{t) (6) 

^ T' 
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Stellar Dynamics around a Massive Black Hole I 7 

is the Hamiltonian, and [ , ](6) is the 6-dim PB defined by| 

df dHors df dHors 


[f , HoTg](Q) — „ • 


(7) 


The Hamiltonian Horg of eqn. 


( 6 ) 

is the sum of an 0(1) term E]^ = (u^/2 — GM,/r^ which is the 
Keplerian orbital energy, and an 0(e) term + r-A,) which comes from the gravity of all the 
stars. Then the CBE can be written as: 

d f 

+ [/) -®k](6) + £[/) + (8) 

Since the Keplerian part is dominant, it is convenient to change the canonical variables from {r, u} 
to new ones, the Delaunay variables in which the Kepler motion of every star appears simplest. 
The Delaunay variables, D = {EL. Ly/, w, q, h}, a re a set of action-angle var iables for the Kepler 

. The three actions 


problem ( Plummer 1960l : Murray &: Dermott 1999 : Binnev Sz Tremaine 


are: I = y/GM,a ; L = ly/l — e^ the magnitude of the angular momentum; and = Lcosi 
the z-component of the angular momentum. The three angles conjugate to them are, respectively: 
w the orbital phase; g the angle to the periapse from the ascending node; and h the longitude 
of the ascending node. When expressed in terms of the D variables, the Keplerian orbital energy 
assumes the simple and degenerate form E\^{I) = —1/2(GM,//)^. Hence for the Kepler problem 
all the Delaunay variables, excepting the orbital phase w, are constant in time; w itself advances 
at the rate 

dEk (GM.)2 


^k(7) - - j3 ■ 

In terms of the D variables, the 6-dim PB between two functions, Xi(T’) and X 2 iD ), is: 

'dxidx2 dxidx2\ , 7dxi dx2 dxi dx2\ , f dxi 9x2 dxi 9x2 


(9) 


[Xl , X2](6) - 


9w 91 91 9w 


+ 


9g 9L dL dg 


+ 


dh dL^ dL. dh 


( 10 ) 

The first PB in eqn. ([8]) is [/, Ek](Q) = Q,\^{df/dw). In the second PB, and A, have to be 
first expressed in terms of the D variables. For this we need to rewrite the phase-space integrals 
of eqns.dS]) and ([3]) in the V variables. The DF is f(V, t) and dr dtr = dV = dl dL dL. dw dq dh . 


Lastly, the position vecto r r{D) = {x,y,z) is given by (|Plummeii 119601 : iMurrav &: Dermottill999l : 
Sambhus fc Sridhar 200oli : 


' X ' 


V 

= 

\ ^ ) 



G„Sh + G^GhS, 


SrS„ 


-SgGh - GiShGg 
-SgSh + GiGhGg 

S,.G„ 


S,Sh \ ( a{Gg -e) \ 


-SiGh 




r\/l — Sg 


0 


( 11 ) 




where S and G are shorthand for sine and cosine of the angles given as subscript. Here a is the 
semi-major axis; rj is the eccentric anomaly, related to the orbital phase through w = {y — e sin y ); 
e = ^l- LyP is the eccentricity; and i is the inclination angle determined by cosi = (L^/L). 
It is also useful to note that r = + y"^ + = a(l — ecosr/). Then eqns. (i2|) and ([3]) take the 

form: 


ip{V,t) = -GM, 


• / tSr’ A.{i) = GM. j 


( 12 ) 


where r = r{T>) and r' = r{V) are given in ean. (|lip . Then the CBE of eqn.([5|l in the D variables 


^ We have used the subscript ‘6’ for the full PB because, in much of this paper, we will use a reduced PB (without subscript) 
— see eqns. lEJ and 12 a — that operates in a 4—dim subspace. 
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8 Sridhar & Touma 


is: 

W + ^•^•](6) = 0- (13) 

This form of the CBE makes apparent two natural time scales in the problem: (i) the Kepler 
orbital period T^ep = 27r/nk , on which the phase w varies, and (ii) the longer secular time scale 
Tsec = £~^Tkep ■ In the next section we orbit-average the CBE over the fast orbital phase w, and 
obtain a reduced CBE describing secular evolution. 


3 ORBIT-AVERAGING VIA THE METHOD OF MULTIPLE SCALES 


Since there are two well-separated time scales (Tkep and Tgec = £' ^Tkep > Tkep) in the problem, 
it is convenient to ana l yze t he CBE of ean. (ll3p using the method of multiple scales (MMS) — see 


e.g 


Bender &: Orszae ( IQTsI l for a general introduction. The MMS proceeds by assuming that all 
relevant quantities — like the DF, potential etc — are functions of two time variables, t and r = et, 
whose variations are treated as independent. This is an artifice whose efficacy is proved when, on 
obtaining an MMS solution, we set t = et to get the desired solution to the full CBE ean. lfT^ : 
thus we write / = /(P, t,r), ip = ip{V,t,T) , A, = A,{t,T). Then, replacing the time derivative 
{d/dt) in the full CBE of equ. lfTS]) by {d/dt + ed/dr) we obtain the MMS-CBE: 


df 


(14) 


df df ^, 

— + + e[/,(^ + r.A.](6) = 0. 

The next step in the MMS is to seek a solution to egn. ([14)1 that describes secular evolution 
plus small fluctuations. Then the dominant part of the DF must be independent of the fast time 
t, as well as the fast orbital phase w. Let us now gather the remaining 5 Delaunay variables and 
call them TZ = {I,L,Lz',g,h} . A point in 5-dim P-space represents a Keplerian ellipse of given 
semi-major axis, eccentricity, inclination, periapse angle and nodal longitude. Henceforth such an 
ellipse will be referred to as a Gaussian Ring or simply Ring. We want to imagine secular evolution 
of A 3> 1 stars as the evolution of A 3> 1 Gaussian Rings over times ~ Tsec • Hence, to 0(1), the 
DF we seek must be a function only of the variables {TZ,t). With these considerations we begin 
with the following ansatz for the DF, as a Fourier expansion in w: 


f{V,t,T) = 


+ 


- y 

27r ^ 


fniTl,t,T)exp[inw] + O(e^), 


(15) 


where T ^ 0 is the 0(1) secular DE, and the 0(e) fluctuation is given by its Fourier-coefhcients 

f-n = fn ■ 

We will derive equations for F and all the fn by substituting the expansion of eon. (1151) in the 
MMS-CBE eqn. ()14l) . and ignoring O(e^) terms. To do this we take the following preparatory steps: 


(a) The quantities, ip{T>, t, r) and A,{t, r), are needed only to 0(1), which implies that we can set 
/ —)> (27r)“^ T(P, r) in eans. lfT^ . 

(b) To 0(1) the self-consistent potential is independent of t, and is given by: 

ifiV, r) = - GM. j F{n', t) dTZ' j> , (16) 

where on the left side we have dropped the dependence on the fast time. Fourier-expanding the 
right side in tc, the potential can be written as: 

ip{'D,T) = 4>(7^,r) + y~)(/?n(7^,r)exp [inw] , (17) 

n^O 
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where 

cI>(7^,r) = 


d 7 ^'F( 7 ^', r) 7 ^'); 


(18a) 


= y d7^'F(7^^r)V'„(7^,7^'), 

with the Ring-Ring interaction potential functions given by 

5.(W,K') = -GAJ.‘ 


J 2tt 2tt \r — r' 


= -CM. ^7“Pi-"”! 


27r 27r 


r — r' 


n / 0 . 


(18b) 

(19a) 

(19b) 


Note that (p-n = V’-n = V^n ) because and 'h are real quantities, 

(c) Using dV = dT^-dtc in the second of eans. lfT^ . we get: 


A, = GM, 


. I F{n\ 


N I drc r 

= O' 


( 20 ) 


because f dwr/r^ = 0 , by the conservation of angular momentum along a Kepler orbit. So we can 
drop the r-A, term in the PB of eon. (1141) . 

(d) The 6 -dim PB of two functions, Xi(D) and X 2 (D), can be expanded as: 

Mxidx2 dxidx2' 


( 1 

[Xl , X 2 j( 6 ) - I - 


dw dl dl dw 
where the PB on the right side (without subscript), 


+ [Xi , X 2 ] , 


( 21 ) 


r 1 ,'dxidx2 dxidx2\ , {dxidx2 dxidx2 
[Xl , X 2 J = 1 - XT-I + 


dg dL dL dg J ' y dh dL^ dL^ dh J ’ 
is the 4-dim PB, whose action is restricted to the 4 dimensional I = constant surfaces in the 5 
dimensional Ring-space. In particular, for two functions Xi{H) and X 2 iH ), the 6 -dim PB equals 
the 4-dim PB: [xi, X 2 ]( 6 ) = [Xi, X 2 ] • 

When eans. (fT5]) and (fT71) are substituted in the MMS-CBE ean. lfT^ . all the 0(1) terms are 
zero. Collecting the 0(e) terms we get: 




(23a) 


df (dF i \ 

+ inOk/^ = + n/0. (23b) 

The second and third terms on the left side of eqn. ()23ap are independent of the fast time t, so 
it is necessary that dfo/dt also be independent of t. Then the general solution is fo{Tl,t,r) = 
Xo{TI,t) + b{JZ,T)t, where Aq and b are arbitrary real functions of TZ and r. However, |/o| is an 
increasing function of t when 6 7 ^ 0 . A cornerstone of the MMS is the requirement that physically 
admissible solutions must not grow on the fast time scale. This implies that 6 = 0, and the MMS 
solution of ean. (|23al) is: 

|^ + [F,4>] = 0, (24a) 


/o(7^,^,r) = Ao(7^,T). (24b) 

Ean. (j24al) is the desired secular CBE for F{TI,t), and is studied in more detail in the next two 
sections. From ean. ()24bh we see that our 0(e) MMS theory only requires /o to independent of t, 
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while letting it be an arbitrary real function of TZ and r. The right side of eQn. (l23bp is independent 
of f, and it is straightforward to write down the general solution: 

1 / dF i \ 

fn{T^,t,T) = + -[F,ipn]j + An(7^, t) exp [-infiki], n/0, (25) 

where the An are arbitrary functions of TZ and r, with A_n = A* . Substituting eans. (j24bp and (1^51) 
in eqn. (IISI), the general solution to the MMS-CBE ean. (jl4p can be written as: 

1 E (dF i \ 


+ -^ 'F “ ^ki)] + O(e^) • 


(26) 


Define two t-independent functions, AF(T>,t) and A{TZ, w,t), by the following Fourier series: 


AF(D, r) = ^ X] ( ] ) exp [inw] , 

n^o ^ ^ 


A{TZ, w,t) = ^ An(7^,r)exp [i 


inw\ 


(27a) 

(27b) 


The Fourier coefficients on the right side of ean. (|27ap are completely determined by F ; they are 
well-defined because n and Dk are both non zero. We do not attempt to prove convergence of the 
Fourier series, and assume that this can be checked in any particular case. In ean. ()27bh the A^ 
are arbitrary functions of TZ and r; they can be chosen such that the Fourier series converges and 
A{TZ, w, t) is a well-defined function of D and r. Having obtained the MMS solution, we go back 
to using a single time variable t, because we want a solution to the original problem of the full 
CBF eqn. (jl3p . In other words, r is no longer regarded as an independent variable, but is simply 
r = et, which is a function of the true time variable t, and the small parameter e. Setting t = et 
in ean. (f26p . we obtain 

f{V,t) = ^F{TZ,et) + ^AF{V,et) + ^A{TZ, w - et) + 0{e^). (28) 

ZTT ZTT ZTT 

which is a solution of the full CBF ean. (ll3p . accurate to 0(e) over times ~ Tgec; this can be verified 
by direct substitution of ean. ([28P in eqn.([T3|). 

The DF of eon. (1280 is the main result of this section. It is the sum of three terms: the first is 
0(1) and proportional to F; the two other terms are 0(e), and can be thought of as the small 
corrections to F. Some salient points to be noted are: 


1. The secular DF F is 0(1), independent of the orbital phase w, and evolves over times Tgec , as 
determined by the CBF ean. (l24ap . 

2. The first 0(e) correction is proportional to AF which is a fluctuating function of the orbital 
phase w with zero mean. It is completely determined by F and, like F, evolves over times Tgec . 

3. The second 0(e) correction is proportional to A which is a function of TZ, {w — Dkt) and et. 
Its functional form is undetermined in our 0{e) MMS theory. However, it is important to note that 
the fast time t occurs only in the form (w — Dkt), which contributes only in a periodic manner, 
with no growth or decay. 

The derivation of the full DF of eqn. (l28p given above makes it natural to think of the orbits 
followed by individual stars, to 0(1) accuracy, as slowly deforming Gaussian Rings. However, it is 
also useful to see this directly from stellar orbital dynamics. This is done in Appendix A using the 
standard apparatus of secular perturbation theory. There it is proved: Every stellar orbit can be 


© 0000 RAS, MNRAS 000 , 000-000 













Stellar Dynamics around a Massive Black Hole I 11 


described to 0(1) accuracy as a Gaussian Ring whose dynamies is governed by the Hamiltonian, 
$(7^, r). Superimposed on this are 0{s) oscillations on the fast time T^ep; so we may think of a 
stellar orbit as a slowly evolving, ‘noisy’ Gaussian Ring. 

MMS analysis is very useful because it deals directly with the DF — bypassing consideration 
of individual stellar orbits — and offers the following: for a given secular evolution described by 
the 0(1) DF F{TZ,t) in 5~dim TZ space, there is a full DF f{'D,t) in the 6~dim V space for which 
the 0(e) noise contributed by all the individual stellar orbits does not act cooperatively to give 
rise to growth on the fast time scale T^ep- Henee the fluetuations remain 0(e) over the long seeular 
time Tgec, and the 0(1) dynamies described by F is well-defined. Thus MMS provides a consistent 
framework for secular dynamics, but it would be well to note that this very consistency has been 
achieved at a certain cost. MMS filters out all 0(1) changes over fast orbital times T^ep, and 
hence the resulting description is insensitive to either fast external perturbations or instabilities. 
It is generally expected that these fast phenomena would have a period of growth followed by 
saturation. Further evolution — including any instability — unfolds over times ~ Tsec- This is the 
regime that is accurately described by the MMS analysis of this Section0 


4 COLLISIONLESS BOLTZMANN EQUATION FOR GAUSSIAN RINGS 


Henceforth we will study the 0(1) slow evolution described by the secular DF F{JZ,t). To this 
order the normalization of the full DF eqn.([T|) implies that 

j F(7^,T)d7^ = 1, (29) 

so we may regard F{TZ,t) as a probability distribution function in 5-dim Ring space. We shall 
refer to F as the Ring distribution function. The time evolution of F is described by the secular 
CBE ean. ()24ap . (dF/dr) + [F, <h] = 0. We recall that is the (scaled) orbit-averaged 

self-gravitational potential of the stars which is determined self-consistently by F itself: 


HH,t) 




(30) 


where 


= -GM, 


dtc du;' 1 
27r 27r \r — r' 


(31) 


is the (scaled) Ring-Ring interaction potential energy. Note that $ is independent of e, so the 
secular evolution of ean. (|24al) has the property of mass-scale invariance: i.e. the behaviour of F 
as a function of TZ and r is independent of £ when e <C 1. This means that, for given MBH mass 
M,, the evolution of the stellar system appears independent of its mass M, so long as M <C M, 
and the time variable used is r = et. This invariance holds when every star is acted upon by the 
Keplerian potentials of the MBH and all the other stars. We now add two corrections to 4>; these 
are the orbit-averaged effects of (a) relativistic corrections and (b) external perturbations. 

(a) Relativistic Corrections to the gravity of the MBH introduce post-Newtonian (PN) cor¬ 
rections, the most important ones being the 1 PN Schwarzschild precession of the apses, and the 
1.5 PN Lense-Thirring precession of the apses and nodes for a spinning MBH. Both effects are 
described by the relativistic (secular) Hamiltonian eR'’®*: 


^ An example of a fast instability is the well-known axisymmetric instability i Toomrelll964 1 applied to a Keplerian disc. All 
but the most unrealistically cold Keplerian discs would be stable by the Toomre criterion. But, as a matter of principle, we 
can indeed consider a Keplerian disc cold enough to be unstable. Then a secular description would make sense only after the 
instability has grown and saturated over times T^ep, leaving the disc in a warmer state that is stable to all fast modes. 
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= hP^\I,L) + 



L) 

II 

1 

'— 1 

CO 1_1 

with 

^ _ 3{GM,)^ M. 

^ M ' 



p 

with 

„ 2{GM,f M, 

(32) 


We have assumed that the spin of the MBH points along the z-axis; 0 ^ cr ^ 1 is the spin 
parameter. Note that the definitions of the constants Bi and ili.s include the factor (M,/M) = . 

(b) External Perturbations on the stars as well as the MBH — due to nuclear density cusps 
and/or slowly moving distant masses — can also be included in the secular description. Let 
r) be the orbit-averaged contribution and be the acceleration of the MBH. 

These will contribute an orbit-averaged tidal potential r) to the secular Hamiltonian for 

a star, where 

+ X(7^)• (33) 


with 

X{n) = f^r{V) (34) 

equal to the orbit-averaged position vector. Note that, similar to the case of relativistic perturba¬ 
tions, has a factor implicit in its definition. 

Putting the pieces together, the general Hamiltonian for a Gaussian Ring is: 

H{n,T) = $(7^,r) + + 4>*''^(7^, r), (35) 


with the Ring equations of motion: 

I = \/GM,a = constant, 

dL _ dH dg_dH dL, _ dH dh _ dH 

dr dg ’ dr dL ^ dr dh ’ dr dL^ 

A secular integral of motion is a function LiJZ, r) that is constant along the orbit 


(36) 


of a 


RingH 


dr 


X(7^(r),r) = 0. 


(37) 


It is a feature of secular dynamics that I = y/GM,a is an integral of motion for quite general and 
time-dependent When the Hamiltonian is time-independent, it provides a second integral 

H{n). 

The equation governing the evolution of the distribution of Gaussian Rings is obtained by 
replacing <1> by 77 in ean. (|24al) : 

^ + [F,H] = 0. (38) 


^ Henceforth we drop the word ‘secular’ when referring to an integral of motion. 


© 0000 RAS, MNRAS 000 , 000-000 












Stellar Dynamics around a Massive Black Hole I 13 


The time evolution is collisionless: using eqns. 
dF _ dF dL dF dg dF , dL^ dF 


dt 


+ 


+ 


dr dL dr dg dr dL 


and (|38p we can verify that 
dh dF 
M~dh 


+ 


dF dH dF 
dr dg dL 

dF 


+ 


dH dF 
dL dg 


dH dF 
dh dL^ 


+ 


dH dF 
dL^ dh 


= ^ + = 0 , 


(39) 

Therefore eon. (1381) describes a system of Gaussian Rings evolving collisionlessly under the combined 
actions of self-gravity, relativistic corrections to the MBH’s Newtonian gravity and slow external 
perturbations. It will be referred to as the Ring CBE. 

The Ring Hamiltonian H of ean. (|35l) is in general not independent of e, because both and 
(^tid proportional to Therefore the property of mass-scale invariance of the Ring CBE is in 
general broken. Another important general property, valid for quite general F{TZ,t) and 
is this: The probability for a Ring to be in (/, I + dl) is a conserved quantity. In other words the 
PDF in I-dim /-space, defined by 


P{I) = j dLdL:,dgdhF{I, L, Lz, g, h, t) , 


(40) 


is independent of r, as can be verified directly using the Ring CBE ean. (f55|l . The physical reason 
for this is that the semi-major axis of every Ring is conserved, so the function P{I) must stay 
frozen as F{TZ,t) evolves over secular times. Therefore secular evolution is driven by gravitational 
torques between Gaussian Rings, with no exchange of Keplerian energies. 


4.1 Collisionless Boltzmann equation for Zero—thickness Flat Discs 

It is useful to study the idealized case when all the Rings are confined to the xy plane. In these 
zero-thickness discs, the angular momentum of every Ring points along Pz, so that a Ring is 
specified by its semi-major axis, angular momentum and apsidal longitude. The Ring phase space 
is now 3-dim: we write TZ = {/, L,g}, where L stands for the angular momentum which can now 
take both positive and negative values {—I ^ L ^ I), and g is the longitude of the periapseH Rings 
with L > 0 will be referred to as prograde and Rings with L < 0 will be referred to as retrograde. 

The general time-dependent Ring DF is F{I, L, g, r). The (scaled) orbit-averaged self-gravitational 
potential, <h(/, L, g(, r), is given in terms of F and d>(TZ,TZ') by eans. dSOl) and (l3T|) . To calculate 
df{TZ,TZ'), we express r = {x,y) in terms of the 2-dim Delaunay variables: 


= a(cos7? — e) cos^f =F ayl — 


e^ sm 7] sm g . 


y = aycos rj — e) sui g 


± aV 1 — 


sm rj cos g . 


(41) 


where the upper(lower) signs in the second terms on the right hand side correspond to pro¬ 
grade (retrograde) orbits. As earlier, rj is related to the mean anomaly hy w = y — esiny. Using 


corresponding expressions for r' = {x',y'), we see that |r —is an explicitly known func¬ 
tion of the Delaunay variables, D and D', and can be averaged over w and w'. Then the secular 
Hamiltonian governing Ring dynamics is: 

H{I,L,g,r) = Hl,L,g,r) + H^^\l,L) + 4>*''^(/, L, ff, r), (42) 


® This is the convention for the symbols (L^g) we will always use when dealing with Ring dynamics that is confined to the xy 
plane. 
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where the relativistic correction, 

L) = -i?i- 


1 


+ 


Sgn(L) 


(43) 


-P\L\ PP 

takes a slightly different form from ean. (|32p . because we have allowed L to take both positive and 
negative values for zero-thickness flat discs. The Ring equations of motion are: 


r clL dH dg dH 

1 = \/GM,a = constant, — = —. 

dr og dr oL 


The DF satisfies the Ring CBE, 

dF dF dH dF 


+ 


dH dF 


dF 


+ [F^H] = 0 , 


(44) 


(45) 


dr dr dg dL dL dg dr 
where the Poisson Bracket is now 2-dim, because it involves only the canonically conjugate pair 
{L,g). As earlier, the probability for a Ring to be in (I, I + dJ) is a conserved quantity: i.e. the 
PDF in 1-dim /-space, defined by 


P{I) = / dLdg F{I,L,g,T), 


(46) 


is independent of r, as can be verified directly using the Ring CBE eon. (14^ . 


5 SECULAR COLLISIONLESS EQUILIBRIA 

Secular collisionless equilibria are described by DFs F that are independent of r; these DFs will 
also be referred to as stationary. From ean. (l38p a stationary DF must satisfy [F,H] = 0. These 
can be constructed using a secular version of Jeans theorem. 


Secular Jeans theorem: Any stationary solution of the Ring CBE ean. (|38l) is a function of TZ 
only through the time-independent integrals of motion of H{TZ), and any function of the time- 
independent integrals of H(TZ) is a stationary so lution of ean.(l38P. 

Proof: This is similar to the one given in Binnev fc Tremaine ( 2008l i: If F’ is a stationary 
solution of ean. (|38l) . then it is also an integral of motion; so the first part of the theorem is proved. 
Conversely, if Xi(7^),... ,XniJZ) be n time-independent integrals, then any F{Xi ,... ,X„) satisfies 
the Ring CBE: 


dr 


because (dXm/dr) = 0 for all m = 1,... 


"F(X,(7^),...,X„(7^)) = Y.dF^. 


m=l 


dim dr 


= 0 , 


(47) 


n. 


□ 


In practice a stronger version of the Secular Jeans theorem is used, by requiring a stationary DF 
to be a function of the isolating integrals of motion. In contrast with non-secular dynamics two 
isolating integrals, I and H{TZ), are always available in secular dynamics. Additional symmetries 
can promote more isolating integrals. Hence a richer variety of stationary DFs can be constructed, 
with greater ease than in non secular dynamics. Below we discuss briefly the simplest cases. 


5.1 Spherical Equilibria 

For a spherically symmetric stellar system around a MBH, the orbit-averaged self-gravitational 
potential 4> is a function only of I and L. An external potential can be included if it is also spherically 
symmetric with respect to the MBH; this causes no acceleration of the MBH, so are 

also functions only of I and L. Therefore the Ring Hamiltonian of eon. (1351) . 

H{fPP) = ck(/,L) + H'^^\I,L,L,) + 4)-t(/,L), (48) 
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is independent of the angle variables. Using eon. (1361) the orbit of a Ring is given by: 
I = constant, L = constant, Lz = constant, 


dr 


n{i,L,Lz) 


dH 


apse precession frequency, 


d/i 

dr 


n{I,L) 


dH 

dLz 


-Bi.5 

(/L)3 ’ 


node precession frequency. 


(49) 


Therefore, every Ring in a spherical system has constant semUmajor axis, eccentricity and in¬ 
clination; its periapse processes at a constant rate U(I, L,L^); its ascending node processes at a 
constant, inclination-independent rate z^(/, L), due solely to 1.5 PN relativistic effects for a spin¬ 
ning MBH. The Ring orbital dynamics is integrable, with the Delaunay actions (I, L,L^) serving 
as three isolating integrals of motion. 

The Secular Jeans theorem implies that a general stationary spherical DF can be a function 
of /, L and Lz ■ We will restrict attention to DFs which have complete spherical symmetry which 
these are of the form T(/, L). R should be noted that the stellar distribution in phase space is given 
explicitly to 0(1), without having to solve an integral equation. Indeed, the functional dependence 
on the six coordinates {r,u} is obtained by setting 

- 1/2 

, L = |r xu| , (50) 

in the DF F[I, L). Then, observationally interesting quantities — such as the radial profiles of the 
surface density, the velocity dispersion tensor and the line-of-sight velocity distribution (LOSVD) 
— can be calculated by direct integration of the DF. It is readily seen that the DFs F{I, L) describe 
systems with zero streaming motions. The velocity dispersion tensor is in general anisotropic. 


I = GM, 


2GM, 9 
- u 


5.2 Axisymmetric Equilibria 

When the stellar system is axisymmetric about the MBH with symmetry axis along z, the orbit- 
averaged self-gravitational potential is a function of (/, L, Lz, g). An external potential can be 
included if it is also axisymmetric with respect to the MBH with the same symmetry axis; this 
causes no acceleration of the MBH, so are also functions only of (/, L, Lz, g). Therefore 

the Ring Hamiltonian of ean. (|35l) simplifies to 

H{I,L,Lz,g) = ^I,L,Lz,g) + H^^\l,L,Lz) + ,L,Lz,g). (51) 

Using ean. (l36l) . the orbit of a Ring is given by 

dg dH dh dH 

I = constant, Lz = constant, — = , — = 7-— . (52) 

dr dL dr dLz 

As a Ring evolves with constant I and Lz , its eccentricity, inclination, apse and node can all vary 
with time. From the Secular Jeans theorem, a general stationary axisymmetric DF must be of the 
form F { I , Lz , H ). Since H = ^ + H '^^^ + is a function of (I, L , Lz , g ), so is the DF. However, 
to express the DF explicitly as a function of the phase space coordinates, we need to determine 
, L , Lz , g ) by solving the self-consistent problem of eans. fpO]) and (f5T]) : 

$(/, L , Lz , g ) = j dl ' dL ' dL '^ dg ' F { l ', L '^, H ') ^ dh ' T(7^, 7^'), (53) 

where H ' = H ( I ', L ' , L '^, g ')^ 

^ The left side of eqn. 1531 1 indicates that $ is independent of /i, as must be the case for any axisymmetric system. We can 
check that the right side of eqn. l53ll is also independent of h\ the nodal longitudes occur in ^(7^,7^') only in the combination 
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The problem simplifies considerably for DFs that are of the form F{I, L^). By replacing I and 
Lz by 

- 1/2 

, Lz = z- (rxit) , (54) 

we see that, jnst like in the spherical case, the DF is known explicitly to 0(1) as a fnnction of 
{r,u}, without having to solve an integral equation. Then observationally interesting quantities 
— such as the sky-maps of the surface density, the velocity dispersion tensor and the line-of-sight 
velocity distribution (LOSVD) — can be calculated by direct integration of the DF. 


I = GM, 


2GM, 


— u 


(a) Zero-thickness Flat Discs: We recall from § 4.1 that the Ring space is now 3-dim; TZ = 
{I,L,g), where —I ^ L ^ I is the angular momentum and g is the apsidal longitude. For an 
axisymmetric disc around a MBH, the (scaled) orbit-averaged self-gravitational potential is a 
function only of I and L. An external potential can be included if it is also axisymmetric with 
respect to the MBH; this causes no acceleration of the MBH, so are also functions 

only of I and L. Therefore the Ring Hamiltonian of eon. (1421) . 

H{I,L) = $(I,L) + + $"^*(I,L), (55) 


is independent of the angle variables. Using ean. (j44p the orbit of a Ring is given by: 

I = constant, L = constant, ^ = Q(I,L) = ——. (56) 

dr oL 

Therefore, every Ring in an axisymmetric system has constant semi-major axis and eccentricity, 
with its apsidal longitude precessing at the constant rate H(/, L). The Ring orbital dynamics is 
integrable, with the Delaunay actions (I, L) serving as two isolating integrals of motion. The Secular 
Jeans theorem implies that a general stationary axisymmetric DF has the form F(I, L). The stellar 
distribution in phase space is given explicitly to 0(1), without having to solve a self-consistent 
problem. Indeed, the functional dependence on the four coordinates {x,y,Ux-,Uy} is obtained by 
setting 


I = GM, 


2GM, 


— u 


- 1/2 


L = xuy - yUx , 


(57) 


in the DF F{I,L). Similar to 3-dim axisymmetric systems, F{I,Lz) discussed above, observa¬ 
tionally interesting quantities — such as the velocity dispersion tensor or the LOSVD — can be 
calculated by direct integration of the DF. Since L can take both positive and negative values, 
a general F (/, L) describes discs with or without rot ation, and generally a nisotr opic velocity dis¬ 
persions. These DFs are relevant to studies such as iTouma Tremaind (j2014l L who solved for 
maximum entropy equilibria of self-gravitating Keplerian disks (in the case where all particles 
have the same semi-major axis). When F(I,L) is an even function of L, local streaming motions 
are absent and the disc is non-rotating with anisotropic velocity dispersions. A special case consists 
of DFs of the form F{I), which are non-rotating and have isotropic velocity dispersions. 


6 PERTURBATIONS AND STABILITY OF SECULAR EQUILIBRIA 

If Fq{TZ) be a stationary DF, the Hamiltonian governing the orbits of Gaussian Rings in the unper¬ 
turbed system is Hq{TZ) = <ho(’^) + L, Lz) + <I>q'^(77), and we must have [Fq ; ^o] = 0. Let 

(h — h'). So, in the integral of $ over h' , when we change the integration variable from h' to {h' — h), we get a quantity that 
is independent of h. 


© 0000 RAS, MNRAS 000 , 000-000 











Stellar Dynamics around a Massive Black Hole I 17 

Fi{TZ,t) be an infinitesimal perturbation about Fq. This gives rise to an infinitesimal perturbation 
of the self-gravitational potential, 




(58) 


Let the external potential also be allowed to change infinitesimally, so it contributes an additional 
r). Then the infinitesimal change in the Ring Hamiltonian is 

tid/ 


= «>i(77,r) + 


(59) 

Substituting F = Fq + Fi and H = Hq + Hi in the Ring CBE eqn. (f38]l . and keeping only the linear 
terms, gives the linearized collisionless Boltzmann equation (LCBE): 

dFi 


dr 


+ [Fi, 77o] + [Fo, iLi] = 0. 


(60) 


Eg ns. ([58)1 — ([60)1 define the linear integral equation for the perturbation Fi(77, r). When = 0, 
the evolution of an infinitesimal perturbation is determined by only its self-gravity $i(77, r) — 
and the unperturbed flows due to Hq{TZ) in which we have allowed for an external contribution 
$0^(77). The unperturbed DE Fo(77) is said to be secularly stable if there is no growing solution 
Fi. Below we discuss two problems: 

(a) Long-lived, warp-like perturbations of spherically symmetric DEs of the form Fo(/, L). 

(b) Axisymmetric, zero-thickness, flat discs with DEs of the form Fq{I,L)- 


6.1 Perturbations of Spherical Equilibria 

Eor a completely spherical unperturbed system the DE is Fq{I,L)j with the unperturbed Ring 
Hamiltonian given by ean. (l48p : 

Ho{I,L,Lz) = ( 61 ) 

Erom eqn.dH]), the apse precession frequency is Q{I, LHz) = dHo/dL; and the inclination- 
independent node precession frequency, zz(/, L) = , is due only to the 1.5 PN relativistic 

correction for a spinning MBH. Then the LCBE of ean. (|60l) can be written as: 

dFo dHi 


dFi BFi dFi 


(62) 


dr dg dh dL dg 

where Hi = contains the effects of the perturbed self-gravity and arbitrary external (slow) 

perturbation. Erom this equation for the linear (secular) perturbation Fi(77, r), we can arrive at 
the following general conclusion: Any growing or decaying behaviour of Fi must necessarily be g- 
dependent, even in the presence of slow external perturbations of arbitrary form. To see this let us 
define the ^-averaged perturbed DE: 


Fi{I,L,Lz,h,T) = j) ^Fi{I,L,Lz,9,h,T), 


(63) 


which can be thought of as the DE for stellar orbits, each of which is smeared into an axisymmetric 
annulus whose inner/outer radii are equal to the peri/apo centre distances. Averaging ean. (l62p 
over g, we get 


m fT - n 

+ zz(i,Fj — 0. 


(64) 


dr ^ ' dh 

The general solution is 

Fi{I,L,Lz,h,r) = Q{I,L,Lz,h- vt) , (65) 

where Q{I,LHz,h) is an arbitrary initial ( 7 -independent perturbation. According to ean. (l 6 ^ such 
an initial condition evolves through precession of the nodes at the rate v{I, F) without either 
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growing or decaying. The nodal precession around the spin axis of the MBH will be negligible in 
non-relativistic stellar systems. Then the perturbations are quasi-static, and could describe long- 
lived warps in spherical nuclear clusters. Since the time evolution — in the linear regime — is 
independent of arbitrary external tidal fields r), the perturbed DF retains memory of its 

origins. 

A linear perturbation must be q-dependent f or it to either grow or decay. We do not inve s tigate 
these, and refer the reader to Tremaine ( 2005 1: IPolvachenko. Polvachenko fc Shukhman (2007) 
which were discussed briefly in the Introduction. Spherical isotropic DFs of the form Fq{I) are 
special because the right hand side of eon. (1621) vanishes. Then the LCBE is: 

TP. a TP. 

( 66 ) 


dFi ^dFi dFi 

+ = 0 . 


dr dg dh 

Therefore an arbitrary infinitesimal perturbation of an isotropic distribution function Fo{I) flows 
collisionlessly along the unperturbed trajectories in 77.-space. The general solution of eqn. dOHI) is: 

Fi{TZ,t) = Q{I, L, Lz, g — h — vt) , forr^O. (67) 

where Q{TZ) = Q{I,L,Lz,g,h) is an arbitrary initial perturbation that is assumed to be given. 
This implies: (a) All isotropic DFs Fo{I) are secularly stable, because the perturbation evidently 
is non growing; (b) The evolution of an initial perturbation is independent of 


6.2 Stability of Axisymmetric Zero—thickness Flat Discs 

We recall from § 4.1 that, for zero-thickness flat discs, the Ring phase space is 3-dim: TZ = 
{I,L,g), where —I ^ L ^ / is the angular momentum, and g is the apsidal longitude. Stationary, 
axisymmetric DFs are of the form Fq{I,L), with the corresponding Hamiltonian of eqn. (l55p : 

Ho{I,L) = ci>o(/,L) - ( 68 ) 


/3|L| PL^ ' u 

From ean. (l56l) we see that, in the unperturbed system, every Ring precesses rigidly at the apse 
precession rate D(/, L) = dHo/dL. A linear perturbation, Fi(77., r), obeys the LCBE: 

dFo dHi 


3F, ^ JF 


8t 


dg 


dL dg ’ 


where 

HiiI,L,g,T) = 4 >,(/,L, 5 ,r) + 4>f(/,L,ff,r), 

cki(/,L, 5 ,r) = fdl'dPdg'^{I,L,g,r,P,g')FiiP,L',g',T) 


(69) 

(70a) 

(70b) 


The linear problem posed by eans. dM]) and (I70bl) is still difficult to solve, because of the integral 
relationship between <ki and Fi in ean. (|70bp . 

Below we derive a stability result, by exploiting the symmetry properties of the Ring-Ring 
interaction potential function d>(TZ,TZ'). We write: 

= -GM. f f. ( 71 ) 

where r = {x,y) and r' = can be written in terms of {I,L,g,g) and {!', L',g', g'), re¬ 

spectively, using ean. (|4T]) : we also recall that dtc = (1 — ecos 7 /)dT/, and dw' = (1 — e'cosT/')d 7 /'. 
The following symmetries of df{TZ,TZ') can be verified either by writing the integral in ean. ()71l) in 
explicit form, or more simply by thinking about the gravitational energy between two Gaussian 
Rings: 

PI: T is a real function which is even in both L and L'. 
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P2: The apsidal longitudes g and g' occur in only in the combination \g — g'\. 

P3: 'h is independent of both g and g' when one of the two Rings is circular (i.e. when L = ±/ or 
L' = ±7' or both). 

P4: is invariant under the interchange of the two Rings. This can be achieved by any of the trans¬ 
formations: {/, L} o {T, , or g g' or both. In explicit form we have: 'b(/, L, T, L', g') = 

'^iI'H',gJ,L,9'), ^{DL,g,r,L',9') = ^{I,L,g',r,L',9), which implies that ^'( 1 , L, 5 ,L', 5 ') = 

^{I',L',g',I,L,9)- 

We now focus attention on disc eigenmodes and hence set = 0 in eqn. ijTOah . so that the 
perturbed Hamiltonian Hi = <hi is due only to the perturbed self-gravity@ Since the coefficients 
in the LCBE of eqn.([69|) are independent of the variables, g and r, it is natural to develop the 
perturbations in a Fourier series. We begin with the Fourier-expansion of the Ring-Ring interaction 
potential function ^{TZ,TZ') in g and g'. Using the fact that 4' depends only on the difference of 
the apsidal longitudes, we write: 

CXD 

'^iI,L,g,l',L',g') = ^ ,L')exv[im{g - g')] . (72) 

m=—oo 

The symmetry properties of 'b, given above in PI to P4, translate to the following properties of 
its Fourier coefficients: 

FI: The Cm are real functions that are even in L and L'. 

F2: The Cm are even in m: i.e. Cm = C-m • 

F3: When either L = ±7 or L' = ±7' or both, then Cm = 0 for all m 7 ^ 0. 

F4: Cm{I,LH'H') = Cm{I',L',DL)- 

The perturbed DF for an eigenmode with azimuthal wavenumber m and eigenfrequency to can be 
written as: 

Fi{I,L, 9 ,r) = Re{Fim(7,L)exp [i(m 5 t - wr)]} . (73) 

where Fi_m = • To compute the corresponding perturbation in the self-gravity, we substitute 

eqns. dZsD and (f72]l in ()70bp . Then the perturbed Hamiltonian is 

Fi( 7 ,L, 5 ,t) = 4 >i( 7 ,L, 5 ,r) = Re{4>im(7,L)exp[i(mc/-wr)]} , 


$1^(7, L) = 27r J dl'dL'Cm{I,L,F,L')Fim{l',L'). 


(74) 


Substituting eans. (|73]) and (|741) in the LCBE eon. ([69]) . we obtain the following linear integral 
eigenvalue problem: 

'dFo' 


[mH(7, L) — oj] Fim = 27rm 


dL 


dl'dH Cm{I,L,l',L')Fim{F,L'), 


(75) 


which must be solved to determine the eigenvalues uj and eigenfunctions Fim{I,L)- 

When m = 0, eqn. (l75j) implies that w = 0 . Therefore all DFs of the form Fq{I, L) are secularly 
stable to axisymmetric perturbations. Moreover these perturbations are also time-independent, 
giving rise to nearby axisymmetric equilibria. It should be noted that this does not imply t hese discs 
are stable to fast axisymmetric instabilities: very cold discs violating the iToomrel (|l964l ) criterion 
will be unstable to axisymmetric modes growing over times ~ Tkep- As we discussed at the end of 
§ 3, MMS analysis is insensitive to these fast modes. Therefore when we consider Keplerian stellar 


® However, in the unperturbed Hamiltonian Hq of ean. ll68l l. the external axisymmetric contribution, L), need not be 

zero. 
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discs, we understand that they must have velocity dispersions that are large enough to make them 
stable to fast modes. 

For non-axisymmetric modes with m 7 ^ 0, the integral eigenvalue problem is, in general, difficult 
to solve. However, we can prove a stability result for a class of stationary DFs. Let Fq{I,L) be a 
strictly monotonic function of L at fixed I. Then we can write 


dFo 

dL 


a 


dFn 


dL 


/ 0 , 


(76) 


where a = +1 when Fq is an increasing function of L, and a = —1 when Fq is a decreasing function 
of L. Let us define a new eigenfunction, Xm{I,L), and a new kernel Am{I, L, I', L') by 


XmiI,L) 

AmiI,L,F,L') 


dFo 

dL 


- 1/2 

Flm{IiL) , 


27r 




dL 


1/2 


C™(/,L,T,L') 


9F' 


dU 


1/2 


(77a) 

(77b) 


where Fq = Fq{I', L'). Then, for m / 0, the eigenvalue eon. (1751) becomes 


H(I,L) 


X^(/,L) 



d/'dL' Am{I,L,l',L')Xm{l',L'). 


(78) 


From F2 and F4 we see that the kernel, Am{I, L, F, L') is real, even in m, and symmetric under 
the interchange {/, L} o {F, L'} . Hence A^ is Hermitian, and all the eigenvalues u are real for 
all m. Therefore we have proved: 


• Stationary, axisymmetric discs with DFs Fq{I,L) are neutrally stable to all secular perturbations 
when dFo/dL is of the same sign (either positive or negative) everywhere in its domain of support, 


I ^ L ^ I and 

-^min 


These secularly stable DFs can have both prograde and retrograde populations of stars because 
-I ^L ^ I. Hence the stability result applies to counter-rotating discs of stars with DFs that are 
strictly monotonic (either increasing or decreasing) functions of L, at fixed /. Depending on the 
sign of dFo/dL, either prograde or retrograde stars dominate at every value of I. These secularly 
stable stellar discs have net rotation and include physically interesting cases, such as a secular 
analogue of the well-known Schwarzschild DF. We may also restate the above stability result as: 


• A necessary condition for Fo{I,L) to he secularly unstable is that dFo/dL must vanish somewhere 
in its domain of support, —I ^ L ^ I and /min ^ / ^ /max • 


Jalali &: Tremaine ( 20121 ) derived the dispersion relationship for slow modes in a collisionless disc 
with a Schwarzschild DF — using the epicyclic approximation and the WKB limit — and showed 
that modes of all m were stable. The stability result presented above is more general, and not 
restricted to a particular DF, or the epi c yclic a pproximation or the WKB limit. Secula r insta - 
bilities have been discussed in Tremaine 


WClic a pproximation or tne Wn.D limit, seculm msta - 

_ ( 2001)1 1: jpolvachenko. Polvachenko fc Shiikhman ( 20071 ). 

Tremaine ( 2005l l considered non-rotating stellar discs described by DFs Fo{I,L) that are 


even 

functi ons of L , and em pty loss cones [Fo/L nonsingular when L —>• 0). Using a variational prin- 
ciple (IG oodma nl 1 198811 he argued that many o f them are likely to be unstable to m = 1 modes. 


Polvachenko. Polvachenko fc Shukhman ( 20071 ) considered mono-energetic counter-rotating discs 


dominated by nearly radial orbits, with DFs equivalent to Fq oc 5{I — Io)A{L) , and found that 
they are prone to loss cone instabilities of all m. 
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7 DISCUSSION 

The main results of this paper are: 

• The full DF in 6-dim phase space, of a secularly evolving, collisionless Keplerian stellar system, 

is given by ean. (|28p . This the sum of an 0(1) secular DF in a reduced 5-dim phase space, and 
small fluctuations that remain of 0(e) over secular times Tsec = ■ 

• The slow time evolution of the secular (or Ring) DF F{TZ,t) in 5-dim 7?.-space is governed 
by the Ring CBE ean. pSI) . which includes the combined secular effects of the gravity of the MBH 
upto 1.5 PN order, the self-gravity of all the stars, and slowly varying external sources of arbitrary 
form. 


The lower dimensional phase space results from the secular conservation of the semi-major axis of 
every stellar orbit. From this fully nonlinear, self-consistent formulation of secular dynamics follows 
a Secular Jeans theorem enabling the construction of Secular Equilibria, and problems concerning 
the Linear Response and Stability of these equilibria. Secular equilibria are easier to construct 
than full stellar dynamical equilibria. For instance, Jeans theorem implies that the exact DF of 
a stationary, spherical, non-rotating system is of the form f{E,L)- In order to determine / as a 
function of the phase space coordinates, {r,u}, it is necessary to solve a self-consistent problem. 
However, the secular DF for the same system, F{I, L), which differs from / only by 0(e), is known 
explicitly as a function of {r,u}. 


_ L inear secular stability of disc-like and spherical configurations has been studied earlier (j Tremaine 

2005 : IPolyachenko, Polyachenko fc Shukhman 2007). Both begin with equilibria and linear pertur¬ 


bations of the full (non-secular) problem. Having formulated an eigenvalue equation for normal 
modes, they then take the slow mode — i.e. a; ~ 0(e) — limit. In contrast we have formulated 
a fully nonlinear secular theory, described by Gaussian Ring dynamics. From this theory followed 
secular equilibria and development of secular linear perturbation theory. Broad conclusions were 
reached about evolutions of small perturbations in spherical non-rotating systems. iTremaind (120051 ) 
has noted that flat systems are less simple than spherical systems due to a difficulty with defining 
inner-products, and the more complicated relation between mass distribution and potential. We 
considered perturbations of DFs describing axisymmetric, zero-thickness flat discs, and derived a 
secular integral eigenvalue problem. Even though the potential is related to the surface density in a 
more complicated manner, its symmetry properties were sufficient to let us prove a stability result 
for a class of physically interesting rotating DFs. 

The fully nonl i near s ecular description of the Ring CBE is needed to go beyond linear theory. 
Touma fc: Sridhar ( 2012l i used a fully nonlinear secular CBE — equivalent to the Ring CBE — to 


discuss the counter-rotating instabil ity, but they did not derive the secular equation from the full 
CBE, as we have done in this paper. Touma, Tremaine &: Kazandiian ( 2009l i performed numerical 
simulations of secular dynamics, wherein each star is replaced by a Gaussian Ring that responds 
to the torques of all the other stars. They found that the counter-rotating instability in unstable 
discs saturates over times ~ Tsec; and the discs settle into uniformly processing equilibria. The 
in-plane instability unfolds into a triaxial one, which couples radial (or eccentricity) and vertical 
(or inclination) degrees-of-freedom, resulting in the dispersal of the lighter of the two counter¬ 
rotating components into a triaxial halo in which the other, more massive, processing, lopsided disc 
is embedded. This process is the secular analogue of “violent relaxation” ( Lvnden-Bell 1967l i and 
forms a central problem in collisionless secular dynamics. The results obtained in this paper are 
used directly in the papers to follow. In Paper H we formulate a theory of the resonant relaxation 
of general Keplerian stellar systems, in terms of the quasi-static evolution (due to gravitational 
collisions between Gaussian Rings) of the Ring DF through a sequence of collisionless equilibria. 
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each of which must be dynamically stable. In Paper III we apply this theory to study the physical 
kinetics of the resonant relaxation of axisymmetric discs. 
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APPENDIX A: SECULAR DYNAMICS OF STELLAR ORBITS 

The full DF, f(T>,t) of eqn. (j28l] . describes the collisionless behaviour of A ^ 1 stars around a 
MBH. The orbit of a test star that belongs to this DF is governed by the Hamiltonian Horg of 
eqn.(j6|l. Horg is the sum of three terms: of these the one containing A, vanishes to 0(e) accuracy 
(see ean J20]l : the Kepler energy is E]^(I) = —1/2 (GM,/!)'^; and the self-gravitational potential is 

-eGM, [ dV = -eGM, [ F{n' ,et) dW I --— 3 -— + 0(e‘^) 

J \r-r'\ J J 2Tr \r - r'\ 

= e^(Tl,et) + e''^^(pn('B,et) exp[mw] + 0(e‘^), (Al) 

n^O 

where we have used eans. ifTOl) and (fT7|l . Hence, to 0(e) accuracy, 

Horg(D,t) = E]^(I) + e^(TZ,et) + e'^^^pn(H,et) exp [inw] , (A2) 

depends on all six Delaunay variables, but has only slow variations with time. Therefore it is natural 
to want to use t = et as the time variable, instead of t. The equations of motion are invariant 
under the rescaling, {t , Horg} —>■ {r = et, Horg = e“^Horg}- We have: 

H'j.g(D,T) = + 4>(77,r) + ^ <^„(77, r) exp [intc] . (A3) 

^ n^O 

The next step makes use of a near-identity canonical transformation from old to new Delaunay 
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variables, V = {!',L',L'^-,w',g',h'}, through the mixed-variable generating function, 

So{l', L', L'/,w, g,h;T) = I'w + L' g + L'^h + eS{l', L', L'/,w, g,h]T) . (A4) 

The old actions {I,L,Lz) and new angles {w\g',h') are given by: 


r r/ dS 

I — I + , 

aw 


r T> 9S 

L — L + , 

dg 


T - T' 


w 


= W + £- 


dS 


9 ' = 9 + £ 


dS 


h’ = h + e 


dS 


(A5) 


5T’ ^ ^ ' “9L'’ 

Hence T> and P', the old and new variables, differ only by 0(e). The new Hamiltonian governing 
the dynamics of T>' is: 

dr 


<g(P',r) = + e- 


dS 

dr 


(A6) 


= (/7„(P, t) exp [inw] + e- 

n^O 

We need to express the right hand side in terms of the new variables V, and keep only the 0(1) 
terms. The first term can be expanded as: 


^k(/) 


= -EJl' + e 


dw 


+ o(£) 


where S^V') = S, L', L'^-, w', g', h'] t) is a function only of the new variables, obtained by re¬ 
placing the old angles {w,g,h) in S by the new ones {w',g',h')\ this is legitimate, because, by 
eqns. (HSl), the old and new variables differ only by 0(e). For the same reason, in the two 0(1) 
terms, we can replace P by Pb Then 

H'org{^',T) = + $(P',t) -k ^¥?n(P',T)exp[inu;'] -k 0(e). (A7) 


n^O 


We now choose S' such that all the tc^-dependent terms on the right hand side cancel. This amounts 
to choosing the mixed variable generating function, 

S{l',L',L'^-,w,g,h;T) = - | ^ P,L(,,ff, fe,r) exp [inw] . (A8) 

S is well-defined, since and n are non-zero. Note that the degeneracy of the Kepler problem 

has prevented the occurrence of small denominators. Then the Hamiltonian for the dynamics of P^ 
is: 


<g(P',r) = + 4>(P',r) + 0(e). 


(A9) 


with the corresponding equations of motion: 


dp 

dr 

= 0(6), — = 

Hk(/') 

e 

' 

■ + gp + Ok). 


dL' 

dr 


dr 

94) 

- au ^ 


dL; 

dr 


dh' 

dr 

94> 

(AlO) 


The dominant behaviour of the orbital phase w' is rapid circulation at the Kepler orbital rate, 
(Hk/^)- Superposed on this are two kinds of modulations; a slow 0(1) variation and a fast 0(e) 
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oscillation. The other five variables, TZ' = {!', L'H'zi tt 'behave like Gaussian Rings with ad¬ 
ditional fast 0(e) oscillations: to 0(1) accuracy, I' is constant, while the other four variables, 
(L', L'zt 9 '^k'), vary over secular time scales with ^(TZ',t) acting as the Hamiltonian for their dy¬ 
namics. This is just the dynamics of TZ' in the orbit-averaged self-gravitational potential. Since, 
by eqn. (jA5p . TZ and TZ' differ only by 0(e), the old TZ variables also have 0(1) slow dynamics 
governed by the Hamiltonian $(7?., r) plus 0(e) fast oscillations. Therefore a stellar orbit in a 
Keplerian stellar system can be thought of as a secularly evolving, “noisy” Gaussian Ring. □ 
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